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We employ a silicon dielectric waveguide to confine and concentrate terahertz pulses, and observe 
that the absorption saturates under strong terahertz fields. By comparing the response between 
lightly-doped and intrinsic silicon waveguides, we confirm the role of hot carriers in this saturable 
absorption. We introduce a nonlinear dynamical model of Drude conductivity that, when incor¬ 
porated into a wave propagation equation, accurately reproduces the observations and elucidates 
the physical mechanisms underlying this nonlinear effect. The results are numerically confirmed by 
Monte Carlo simulations of the Boltzmann transport equation, coupled with split-step nonlinear 
wave propagation. 


Among semiconductors, silicon is not only the most 
prevalent material in electronics, but it is also one of 
the most favorable dielectric materials for terahertz ap¬ 
plications. Intrinsic silicon is transparent at wavelengths 
longer than 1100 nm, and has exceptionally low loss in 
the far infrared [T]. While the nonlinear properties and 
applications of silicon are well established in the near- 
infrared and mid-infrared regime [2], there have been very 
few observations of nonlinear propagation in the tera¬ 
hertz regime. 

The terahertz photon energy (4.1 meV at 1 THz) is too 
small to produce new carriers in silicon through a 1- or 
2-photon absorption, and hence the linear and nonlinear 
properties are caused by acceleration or heating of the 
existing electron (or hole) population. The traditional 
Drude model of conductivity that is commonly used to 
describe free carrier absorption and dispersion in silicon 
in the terahertz regime fails to explain nonlinear wave 
propagation effects. 

In 2010, Hebling et al. and Kaur et al. indepen¬ 
dently observed THz field induced absorption bleach¬ 
ing in n-doped bulk silicon, using terahertz pump-probe 
measurements [3] and z-scan measurements [4]. They sug¬ 
gested that the effect might be explained by scattering of 
electrons into a higher energy (L) valley within the con¬ 
duction band. Terahertz induced nonlinear effects have 
also been observed in a variety of other bulk semiconduc¬ 
tors, including Ge[3], GaAs[5], GaP[6] and InSb[7], and 
numerous hot carrier effects have been offered as explana¬ 
tions, including intervalley scattering, band nonparabol- 
icity, and impact ionization. In most cases, the obser¬ 
vations were carried out using wafers or windows with 
optical thickness of only a few terahertz wavelengths. In 
such thin samples, the cumulative nonlinearity is neces¬ 
sarily quite small, and it is difficult to separate propaga¬ 
tion effects from interface effects such as small changes 
in reflectivity, or spatial effects such as self-focusing and 
diffraction. 

To overcome these limitations, we couple picosecond 
terahertz pulses into a 2 cm long silicon dielectric ridge 
waveguide. The waveguide greatly enhances the field 


concentration and nonlinear propagation length, thereby 
ensuring that the measured effect represents a true non¬ 
linear wave interaction accumulated over hundreds of 
terahertz wavelengths, and also allows for interplay be¬ 
tween the linear mode propagation and nonlinearity. The 
waveguide configuration also eliminates spatial nonlinear 
effects like self-focusing, enabling unambiguous measure¬ 
ment of the temporal nonlinear behavior. We observe 
over a two-fold increase in the power transmission ratio 
at high powers relative to low powers, depending on the 
carrier concentration, and we present a new physical and 
numerical model that explains the observed behavior. 

The silicon ridge waveguides were fabricated from 
400 jim thick, double-side polished (DSP), (100) silicon 
wafers. In order to better assess the role of carriers, we 
used two types of silicon: lightly p-doped wafers with 
a nominal resistivity of of 150-350 Q-cm and float-zone 
semi-insulating wafers with a resistivity of 10 kf^-cm. A 
1 fim sacrificial layer of Si02 was deposited by CVD on 
the wafers, and patterned using contact photolithography 
and reactive-ion etching to produce a 300 jam wide oxide 
hard-mask for subsequent etching of the waveguides. The 
waveguides were etched to a depth of 100 fim using pulsed 
deep reactive ion etching (Bosch process), after which 
the remaining photoresist and oxide hard mask were re¬ 
moved. Fig. [^a) shows a cross-sectional micrograph of 
the completed ridge waveguide, and Fig. Bb) shows the 
corresponding fundamental TE eigenmode of the waveg- 



FIG. 1. (a) Cross-sectional micrograph of fabricated silicon 
ridge waveguide and (b) calculated TE eigenmode at 0.5 THz. 
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FIG. 2. Experimental setup used to measure the THz nonlin¬ 
ear transmission through the silicon waveguide. 


FIG. 3. Normalized power transmission for semi-insulating 
(circle) and doped (square) waveguides, and corresponding 
calculated (dashed lines) pulse energy transmission. 


uide, calculated at 0.5 THz. The transverse waveguide 
dimensions were chosen to ensure single-mode operation 
over the frequency range of interest. The waveguides 
were cut to a length of 2 cm using a dicing saw. 

Fig-i illustrates the experimental setup used to char¬ 
acterize the THz nonlinear response. An amplified Ti- 
sapphire laser system produces 40 fs, 1 kHz repetition 
rate pulses at 800 nm center wavelength. The optical 
pulses are split (80:20) into pump and probe beams that 
are used for terahertz generation and detection, respec¬ 
tively. The pump pulse impinges on a grating (2000 
lines/mm), producing a —1 order diffracted beam that 
has a tilted pulse front [8]. The tilted pulse was de- 
magnified by a factor of 2x using a 60 mm focal length 
lens into a LiNbOs prism. A A/2 waveplate rotates the 
optical beam polarization from horizontal to vertical di¬ 
rection to align with the optical axis of the LiNbOs. 
The power of the THz output beam was adjusted us¬ 
ing a pair of wire-grid polarizers, and focused using a 
polymethylpentene (TPX) lens onto the input waveguide 
facet. The THz beam was linearly polarized in the (011) 
crystallographic direction of the silicon waveguide. Using 
the experimentally measured energy, pulse duration, and 
focused spot size of the terahertz beam, the peak elec¬ 
tric field at the focus before inserting the waveguide was 
estimated to be 200 kV/cm[8]. 

The THz pulses impinging on and emerging from the 
waveguides were measured using both a pyroelectric de¬ 
tector and electrooptic sampling. In the latter case, we 
used a 1 mm thick (110) ZnTe crystal that was coated 
with an 800 nm dielectric mirror front face, and antire¬ 
flection coating on the rear face, which allows the probe 
beam to be introduced in a reflection geometry [9], as 
shown in Fig.|^ The ZnTe electrooptic crystal was placed 
in contact with the output facet of the waveguide, to al¬ 
low for near-held optical sampling of the mode emerging 
from the waveguide. 

To measure the nonlinear transmission through the 


waveguide, we used the Fourier transform to calculate the 
spectrum of the emerging waveform, and integrated the 
intensity spectrum to obtain a measure of the transmit¬ 
ted power. For the range of powers considered, the non¬ 
linearity of the electrooptic detection process was con- 
hrmed to be negligible in comparison to the absorption 
saturation in the silicon waveguide. 

Fig-i shows the normalized transmission ratio as a 
function of the input pulse energy and peak held for the 
two waveguides considered here. The semi-insulating sil¬ 
icon waveguide shows a small, but clearly measurable 
increase of 5% in transmission as the pulse energy is in¬ 
creased from 0 to 75 nJ. The p-type silicon waveguide, by 
contrast, shows a more than 2-fold increase in transmis¬ 
sion at higher huence. The dashed lines plot the numer¬ 
ically calculated result (to be explained below), which 
shows that at sufficiently high pulse power, the power 
transmission ratio saturates at a level close to unity. The 
fact that the saturable absorption is much stronger in 
doped silicon clearly demonstrates the role of free carri¬ 
ers in the nonlinear response. 

A complete model of absorption in silicon waveguides 
must account for not only the field-dependent nonlinear 
carrier dynamics, but also the linear dispersion, which 
diminishes the peak field of the signal. The terahertz 
nonlinear wave propagation can be described by a sim¬ 
plified one-dimensional wave equation. 
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where J is the current density (which is non-linearly re¬ 
lated to E) and P is the linear polarization of the mate¬ 
rial, which is linearly related to the electric field in the 
frequency domain by: 

p{z,uj) = eo N(w) - l] E{z,uj ), 


( 2 ) 
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where n{uj) is the effective the refractive index of the 
waveguide. 

If the current is neglected, the forward traveling solu¬ 
tion to 0 in the frequency domain is 


E{Az, uj) = E{0^ uj) exp 


/ N A ' 

i—niuj)Az 
. c 


( 3 ) 


where the refractive index n{uj) incorporates material and 
modal dispersion of the waveguide. 

Conversely, if the dispersion is neglected but the cur¬ 
rent term is retained, then the wave equation can be writ¬ 
ten as: 


d^E 1 d^E _ dJ 

where v = c/n{Co) represents the average velocity of the 
terahertz pulse, evaluated at the center frequency of the 
spectrum. 

Because the electric field travels in the -\-z direction 
with an average velocity of i;, we assume that the result¬ 
ing current density can be likewise cast as a function of 
a single argument, J{t — z/v)^ in which case 0 can be 
integrated to find the field at the end of one step A 2 :: 


E{Az^ t) = E'(0, t — Az/v) 

t—Azjv 

+ ^ / [J{t')-J{t-^z/v)]dt' (5) 

The second term in ^ represents a perturbation AE in 
the electric field caused by the current J. The split-step 
numerical method replaces this accumulated nonlinearity 
by an equivalent lumped effect at 2 : = 0, which is found 
by advancing ^ by the propagation time Azjv^ 

t 

J [Jit')-J{t)]dt' (6) 

t—2Azlv 


The nonlinear wave propagation is numerically sim¬ 
ulated by dividing the total propagation distance into 
steps of size Az, computing the linear propagation for 
each increment in the Fourier domain using Q, and in¬ 
corporating the nonlinearity as lumped in the time do¬ 
main using 0. 

The nonlinear relationship between the electric field 
E{t) and current density J{t) can be described using the 
balance equations obtained from the Boltzmann trans¬ 
port equations. In the spatially homogeneous limit, the 
momentum balance equation is m¬ 


dv 

dt 




m* 


( 7 ) 


where v represents the carrier velocity, which is directly 
proportional to the current density through J = Nqv^ 


and Tmis) is the momentum relaxation rate, which we 
take to be a function of the energy, 5 . 

The energy balance equation is 

dr 

-+T,e = qEv, (8) 

where 5 is the carrier energy relative to thermal equilib¬ 
rium, and Tg is the energy relaxation rate. The momen¬ 
tum and energy scattering rates are, in general, energy 
dependent, which couples these two equations. We adopt 
the simple, and widely used model where the energy re¬ 
laxation rate T^ is taken to be constant, while the mo¬ 
mentum relaxation rate increases linearly with the carrier 
energy [11]: 


r„(e) = To + (9) 

For sufficiently small carrier energy, the second term in 
0 may be neglected, in which case 0 can be solved 
directly to give the familiar linear Drude relationship be¬ 
tween V and E, in the frequency domain, 

viu;)= ( 10 ) 

1 — 100/ L 0 

where /i = q/mTo is the low-field mobility. However, for 
sufficiently high fields, Eqs. 0-0 predict well-known 
nonlinear transport phenomena including the saturation 
of carrier velocity at v = with increasing DC field 
strength. The electron and hole saturation velocities in 
silicon are approximately i;sat ^ 10^ cm/s, and the corre¬ 
sponding critical electric field strength above which sat¬ 
uration effects become important is Ecr = '^sat/M ^ 7 
kV/cm (for electrons), 16 kV/cm (for holes) - conditions 
that are readily achieved for the terahertz pulses used in 
these experiments. 

We used the split-step numerical method described 
above, together with the nonlinear Drude relations de¬ 
scribed in 0-0 to calculate the power-dependent trans¬ 
mission as a function of input power for the two waveg¬ 
uides under consideration. For the p-doped silicon sam¬ 
ple, we assumed a carrier concentration of N = 8.5 x 10^^ 
cm“^, a low-field hole mobility of /i = 470 cm^/(V-s), 
a hole effective mass of m* = 0.36mo, and a satura¬ 
tion velocity of 0.75 x 10^ cm/s. For the high-resistivity 
fioat-zone silicon, we estimated a residual electron con¬ 
centration of N = 5 X 10^^ cm“^ and a low-field elec¬ 
tron mobility of /i = 1,416 cm^/(V's), effective mass 
of m* = 0.26mo, and an electron saturation velocity of 
10^ cm/s. In both cases the energy relaxation rate was 
taken to be l/F^ = 0.2 ps. We used accepted physical 
parameters from the literature, and the only adjustable 
parameter in the calculation was the carrier concentra¬ 
tion A, which was chosen to both match the resistivity 
range of the wafers and to also agree with the observed 
absorption in the low-field limit. The calculations were 
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(a) (b) 

FIG. 4. Evolution of the temporal profile of terahertz pulse 
along the waveguide obtained from (a)conventional Drude 
model and (b)nonlinear split-step simulation 

performed using 100 /im steps and a temporal window 
of 80 ps divided into 4000 steps. For the numerical cal¬ 
culations, the input THz waveform was taken to be of 
the form E{t) = Eq cos{at)e~^^ where the constants a 
and b were chosen to best match the actual measured 
input waveform. The dashed lines in Fig. show the 
calculated pulse energy transmission as a function of the 
input power and peak field, and agree well with the ex¬ 
perimental measurements. 

Fig. 11 shows a numerical simulation of how the tera¬ 
hertz pulse evolves in time as it traverses the 2 cm long 
waveguide. The left portion was calculated using the con¬ 
ventional Drude model while the right portion includes 
the nonlinear split-step model discussed here, assuming a 
peak-peak input field of 100 kV/cm, clearly showing the 
enhanced field transmission. 

The balance equations 0-® provide a simple and ef¬ 
ficient model of the nonlinear transport in silicon, but 
there are alternative models used to explain the energy- 
dependent relaxation rates. The most accurate and 
widely accepted approach is to use the Monte Carlo 
method to directly simulate the Boltzmann transport 
equations in the time domain [12]. To better resolve the 
physical origins of the nonlinearity, we used the same 
split-step Fourier method to compute the nonlinear prop¬ 
agation, but instead of 0-®, the current density at 
each step was estimated using time-dependent Monte 
Carlo simulations of an ensemble of 10,000 carriers. This 
method is far more computationally intensive, and we 
therefore divided the waveguide into only 20 steps and 
simulated propagation for an input pulse with peak-peak 
field of 100 kV/cm. The same enhancement of transmis¬ 
sion is observed (Fig. [^. 

The Monte Carlo calculations incorporate several 
physical effects that contribute to the observed response, 
including band non-parabolicity. Coulomb scattering, in¬ 
travalley acoustic phonon scattering, and equivalent in¬ 
tervalley optical phonon scattering. Of these, simula¬ 
tions revealed that intravalley and equivalent intervalley 




FIG. 5. (a) Transmitted terahertz waveform, calculated using 
Monte Garlo simulation of carrier dynamics together with the 
split-step Fourier method. The linear (green) output curve 
was calculated using the conventional linear Drude model 
and waveguide dispersion, and shows lower transmission. In¬ 
set: the simulated input pulse, with a peak-peak field of 100 
kV/cm. The Supplemental Material at [URL will be inserted 
by publisher] provides an animation showing how the field 
and nonlinearity evolve with distance, (b) Transmitted power 
spectrum, calculated with (blue) and without (green) nonlin¬ 
earity. (c) Experimentally measured transmitted terahertz 
waveform for 75 nJ (blue) and 0.75 nJ (green) incident pulse 
energy. The green field was scaled by 10 x to account for the 
100 X lower energy, (d) Experimentally measured power spec¬ 
tra for 75 nJ (blue) and 0.75 nJ (green) incident pulse energy. 
The green curve was scaled by 100 x to account for the lower 
energy. 

phonon scattering were found to be the dominant factors 
that contribute to the nonlinearity in the simulated re¬ 
sponse. Notably, higher energy L-X intervalley scattering 
does not play a significant role, as had been previously 
suggested. 

Fig.[^a)-(b) show the calculated output waveform and 
spectra, obtained using a combination of the Monte Carlo 
simulation with split step Fourier method, for the high¬ 
est input field (100 kV/cm) that was considered in the p- 
doped waveguide. For comparison we also show the field 
obtained from the conventional (linear) Drude model, 
which would predict a higher carrier velocity and lower 
output field. The Supplemental Material at [URL will 
be inserted by publisher] provides an animation showing 
how the nonlinearity and dispersion accumulate as the 
pulse traverses waveguide. 

Fig. i c)-(d) show the corresponding experimental 
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measurements of the output terahertz waveforms and 
spectra, which show a similar increase in relative trans¬ 
mission at high fields. The experimental spectra show 
additional loss that is attributed to strong water absorp¬ 
tion at 0.55, 0.75 and 1.1 THz that is absent from the 
simulations. To assess the role of nonlinearity, we atten¬ 
uated the input power by a factor of 100 x and repeated 
the measurement of the output waveform. The green lin¬ 
ear curve shown in Fig.[^c)-(d) was then scaled by a fac¬ 
tor of 10 or 100 to provide a direct comparison with the 
field and power (respectively) measured at higher power. 
These measurements clearly show that the relative trans¬ 
mission is significantly higher for strong terahertz pulses, 
and the degree of absorption saturation is comparable to 
that shown in Fig. [^b). 

In conclusion, we experimentally explore the phe¬ 
nomenon of absorption saturation in silicon dielectric 
waveguides at terahertz frequencies. The field-induced 
transparency and associated carrier velocity saturation 
is shown to be dynamical effect that cannot be ade¬ 
quately explained by a modified effective mobility or 
Drude model. We present a simple, nonlinear Drude 
model that explains the observations, and we confirm the 
model using rigorous Monte Carlo simulations. Further, 
we introduce a numerical split-step method that mod¬ 
els the interplay of nonlinearity and dispersion in the 
wave propagation. These results could have important 
consequences in future high-power terahertz guided-wave 
nonlinear devices, such as terahertz frequency converters, 
parametric oscillators, mixers, and modulators. 
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